

/*
 * Class:        BrownianMotionBridge
 * Description:  
 * Environment:  Java
 * Software:     SSJ 
 * Copyright (C) 2001  Pierre L'Ecuyer and Universite de Montreal
 * Organization: DIRO, Universite de Montreal
 * @author       
 * @since

 * SSJ is free software: you can redistribute it and/or modify it under
 * the terms of the GNU General Public License (GPL) as published by the
 * Free Software Foundation, either version 3 of the License, or
 * any later version.

 * SSJ is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.

 * A copy of the GNU General Public License is available at
   <a href="http://www.gnu.org/licenses">GPL licence site</a>.
 */

package umontreal.iro.lecuyer.stochprocess;
import umontreal.iro.lecuyer.rng.*;
import umontreal.iro.lecuyer.probdist.*;
import umontreal.iro.lecuyer.randvar.*;



/**
 * Represents a Brownian motion process 
 * <SPAN CLASS="MATH">{<I>X</I>(<I>t</I>) : <I>t</I>&nbsp;&gt;=&nbsp;0}</SPAN>
 * sampled using the <SPAN  CLASS="textit">bridge sampling</SPAN> technique
 * (see for example).
 * This technique generates first the value <SPAN CLASS="MATH"><I>X</I>(<I>t</I><SUB>d</SUB>)</SPAN> at the last observation time,
 * then the value at time <SPAN CLASS="MATH"><I>t</I><SUB>d/2</SUB></SPAN> (or the nearest integer),
 * then the values at time <SPAN CLASS="MATH"><I>t</I><SUB>d/4</SUB></SPAN> and at time <SPAN CLASS="MATH"><I>t</I><SUB>3d/4</SUB></SPAN>
 * (or the nearest integers), and so on.
 * If the process has already been sampled at times <SPAN CLASS="MATH"><I>t</I><SUB>i</SUB> &lt; <I>t</I><SUB>k</SUB></SPAN> but not
 * in between, the next sampling point in that interval will be
 * <SPAN CLASS="MATH"><I>t</I><SUB>j</SUB></SPAN> where 
 * <SPAN CLASS="MATH"><I>j</I> = floor((<I>i</I> + <I>k</I>)/2)</SPAN>.
 * For example, if the sampling times used are
 * {
 * <SPAN CLASS="MATH"><I>t</I><SUB>0</SUB>, <I>t</I><SUB>1</SUB>, <I>t</I><SUB>2</SUB>, <I>t</I><SUB>3</SUB>, <I>t</I><SUB>4</SUB>, <I>t</I><SUB>5</SUB></SPAN>},
 * then the observations are generated in the following order:
 * <SPAN CLASS="MATH"><I>X</I>(<I>t</I><SUB>5</SUB>)</SPAN>, <SPAN CLASS="MATH"><I>X</I>(<I>t</I><SUB>2</SUB>)</SPAN>, <SPAN CLASS="MATH"><I>X</I>(<I>t</I><SUB>1</SUB>)</SPAN>, <SPAN CLASS="MATH"><I>X</I>(<I>t</I><SUB>3</SUB>)</SPAN>, <SPAN CLASS="MATH"><I>X</I>(<I>t</I><SUB>4</SUB>)</SPAN>.
 * 
 * <P>
 * <SPAN  CLASS="textit">Warning</SPAN>:
 * Both the <TT>generatePath</TT> and the <TT>nextObservation</TT> methods from
 * {@link umontreal.iro.lecuyer.stochprocess.BrownianMotion BrownianMotion} are
 * modified to use the bridge method.
 * 
 * In the case of <TT>nextObservation</TT>, the user should understand
 * that the observations returned are <SPAN  CLASS="textit">not</SPAN> ordered chronologically.
 * However they will be once an entire path is generated and the observations
 * are read from the internal array (referenced by the <TT>getPath</TT> method)
 * that contains them.
 * 
 * <P>
 * The method <TT>nextObservation(double nextTime)</TT> differs from that of
 * the class
 * {@link umontreal.iro.lecuyer.stochprocess.BrownianMotion BrownianMotion}
 *  in that <TT>nextTime</TT> represents
 * the next observation time <SPAN  CLASS="textit">of the Brownian bridge</SPAN>.
 * However, the <SPAN CLASS="MATH"><I>t</I><SUB>i</SUB></SPAN> supplied must still be non-decreasing with <SPAN CLASS="MATH"><I>i</I></SPAN>.
 * 
 * <P>
 * Note also that, if the path is not entirely generated before being read
 * from this array, there will be ``pollution'' from the previous path generated,
 * and the observations will not represent a sample path of this process.
 * 
 */
public class BrownianMotionBridge extends BrownianMotion  {
    protected int          bridgeCounter = -1; // Before 1st observ

    // For precomputations for B Bridge
    protected double[]     wMuDt,
                           wSqrtDt;
    protected int[]        wIndexList,
                           ptIndex;



   /**
    * Constructs a new <TT>BrownianMotionBridge</TT> with
    * parameters 
    * <SPAN CLASS="MATH"><I>&#956;</I> = <texttt>mu</texttt></SPAN>, 
    * <SPAN CLASS="MATH"><I>&#963;</I> = <texttt>sigma</texttt></SPAN> and initial value
    * 
    * <SPAN CLASS="MATH"><I>X</I>(<I>t</I><SUB>0</SUB>) = <texttt>x0</texttt></SPAN>.
    * The normal variates will be
    * generated by inversion using the
    * {@link umontreal.iro.lecuyer.rng.RandomStream RandomStream} <TT>stream</TT>.
    * 
    */
   public BrownianMotionBridge (double x0, double mu, double sigma,
                                RandomStream stream)  {
        super (x0, mu, sigma, stream);
    }


   /**
    * Constructs a new <TT>BrownianMotionBridge</TT> with
    * parameters 
    * <SPAN CLASS="MATH"><I>&#956;</I> = <texttt>mu</texttt></SPAN>, 
    * <SPAN CLASS="MATH"><I>&#963;</I> = <texttt>sigma</texttt></SPAN> and initial value
    * 
    * <SPAN CLASS="MATH"><I>X</I>(<I>t</I><SUB>0</SUB>) = <texttt>x0</texttt></SPAN>.
    * The normal variates will be
    * generated by the
    * {@link umontreal.iro.lecuyer.randvar.NormalGen NormalGen} <TT>gen</TT>.
    * 
    * 
    */
   public BrownianMotionBridge (double x0, double mu, double sigma,
                                NormalGen gen)  {
        super (x0, mu, sigma, gen);
    }


   public double nextObservation() {
        double x;
        if (bridgeCounter == -1) {
            x = x0 + mu*(t[d]-t[0]) + wSqrtDt[0] * gen.nextDouble ();
            bridgeCounter = 0;
            observationIndex = d;
        } else {
           int j = bridgeCounter*3;
           int oldIndexL = wIndexList[j];
           int newIndex  = wIndexList[j + 1];
           int oldIndexR = wIndexList[j + 2];

           x = path[oldIndexL] +
             (path[oldIndexR] - path[oldIndexL])
             * wMuDt[newIndex] + wSqrtDt[newIndex] * gen.nextDouble ();

           bridgeCounter++;
           observationIndex = newIndex;
        }
        observationCounter = bridgeCounter + 1;
        path[observationIndex] = x;
        return x;
    }

   public double nextObservation (double nextTime) {
        double x;
        if (bridgeCounter == -1) {
            t[d] = nextTime;

            wMuDt[0]   = 0.0;  // The end point of the Wiener process
                               //  w/ Brownian bridge has expectation = 0
            wSqrtDt[0] = sigma * Math.sqrt(t[d] - t[0]);
                               // = sigma*sqrt(Dt) of end point

            x = x0 + mu*(t[d]-t[0]) + wSqrtDt[0] * gen.nextDouble ();

            bridgeCounter = 0;
            observationIndex = d;
        } else {
            int j = bridgeCounter*3;
            int oldIndexL = wIndexList[j];
            int newIndex  = wIndexList[j + 1];
            int oldIndexR = wIndexList[j + 2];

            t[newIndex] = nextTime;

            double dtRL = t[oldIndexR] - t[oldIndexL];
            if (dtRL != 0.0) {
                wMuDt[newIndex] = (t[newIndex]-t[oldIndexL]) / dtRL;
            } else {
                wMuDt[newIndex] = 0.0;
            }
            wSqrtDt[newIndex] = sigma * Math.sqrt (
               wMuDt[newIndex] * (t[oldIndexR] - t[newIndex]));

            x = path[oldIndexL] +
              (path[oldIndexR] - path[oldIndexL])
              * wMuDt[newIndex] + wSqrtDt[newIndex] * gen.nextDouble ();

            bridgeCounter++;
            observationIndex = newIndex;
        }
        observationCounter = bridgeCounter + 1;
        path[observationIndex] = x;
        return x;
    }

   public double[] generatePath() {
        // Generation of Brownian bridge process
        int oldIndexL, oldIndexR, newIndex;
        path[d] = x0 + mu*(t[d]-t[0]) + wSqrtDt[0] * gen.nextDouble ();

        for (int j = 0; j < 3*(d-1); j+=3) {
           oldIndexL   = wIndexList[j];
           newIndex    = wIndexList[j + 1];
           oldIndexR   = wIndexList[j + 2];

           path[newIndex] = path[oldIndexL] +
             (path[oldIndexR] - path[oldIndexL])
             * wMuDt[newIndex] + wSqrtDt[newIndex] * gen.nextDouble ();
        }

        //  resetStartProcess();
        observationIndex   = d;
        observationCounter = d;
        return path;
    }

    public double[] generatePath (double[] uniform01){
        int oldIndexL, oldIndexR, newIndex;
        path[d] = x0 + mu*(t[d]-t[0]) + wSqrtDt[0] * NormalDist.inverseF01(uniform01[0]);

        for (int j = 0; j < 3*(d-1); j+=3) {
           oldIndexL   = wIndexList[j];
           newIndex    = wIndexList[j + 1];
           oldIndexR   = wIndexList[j + 2];

           path[newIndex] = path[oldIndexL] +
             (path[oldIndexR] - path[oldIndexL])
             * wMuDt[newIndex] + wSqrtDt[newIndex] * NormalDist.inverseF01(uniform01[1 + j/3]);
        }

        //  resetStartProcess();
        observationIndex   = d;
        observationCounter = d;
        return path;
    }

    public void resetStartProcess() {
        observationIndex   = 0;
        observationCounter = 0;
        bridgeCounter = -1;
    }

    protected void init() {
      super.init();

      /* For Brownian Bridge */

      // Quantities for Brownian Bridge process
      wMuDt = new double[d + 1];
      wSqrtDt = new double[d + 1];
      wIndexList = new int[3 * (d)];
      ptIndex = new int[d + 1];
      double tem = 0;

      int indexCounter = 0;
      int newIndex, oldLeft, oldRight;

      ptIndex[0] = 0;
      ptIndex[1] = d;

      wMuDt[0] = 0.0;  // The end point of the Wiener process
      //  w/ Brownian bridge has expectation = 0
      if (t[d] < t[0])
         throw new IllegalStateException("   t[d] < t[0]");
      wSqrtDt[0] = sigma * Math.sqrt(t[d] - t[0]);
      // = sigma*sqrt(Dt) of end point

      for (int powOfTwo = 1; powOfTwo <= d / 2; powOfTwo *= 2) {
         /* Make room in the indexing array "ptIndex" */
         for (int j = powOfTwo; j >= 1; j--) {
            ptIndex[2*j] = ptIndex[j];
         }

         /* Insert new indices and Calculate constants */
         for (int j = 1; j <= powOfTwo; j++) {
            oldLeft = 2 * j - 2;
            oldRight = 2 * j;
            newIndex = (int) (0.5 * (ptIndex[oldLeft] + ptIndex[oldRight]));

            wMuDt[newIndex] = (t[newIndex] - t[ptIndex[oldLeft]]) /
                              (t[ptIndex[oldRight]] - t[ptIndex[oldLeft]]);
            tem = (t[newIndex] - t[ptIndex[oldLeft]]) *
                  (t[ptIndex[oldRight]] - t[newIndex])
                  / (t[ptIndex[oldRight]] - t[ptIndex[oldLeft]]);

            // Test for NaN (z != z); 0/0 gives a NaN
            if (tem < 0 || tem != tem) {
               System.out.printf ("t[newIndex] - t[ptIndex[oldLeft]] = %g%n", t[newIndex] - t[ptIndex[oldLeft]]);
               System.out.printf ("t[ptIndex[oldRight]] - t[newIndex] = %g%n", t[ptIndex[oldRight]] - t[newIndex]);
               System.out.printf ("t[ptIndex[oldRight]] - t[ptIndex[oldLeft]] = %g%n", t[ptIndex[oldRight]] - t[ptIndex[oldLeft]]);
               System.out.printf ("t[ptIndex[oldRight]] = %g%n", t[ptIndex[oldRight]]);
               System.out.printf ("t[ptIndex[oldLeft]] = %g%n", t[ptIndex[oldLeft]]);
               throw new IllegalStateException("   tem < 0 or NaN");
            }
            wSqrtDt[newIndex] = sigma * Math.sqrt (tem);

            ptIndex[oldLeft + 1] = newIndex;
            wIndexList[indexCounter] = ptIndex[oldLeft];
            wIndexList[indexCounter + 1] = newIndex;
            wIndexList[indexCounter + 2] = ptIndex[oldRight];

            indexCounter += 3;
         }
      }
      /* Check if there are holes remaining and fill them */
      for (int k = 1; k < d; k++) {
         if (ptIndex[k - 1] + 1 < ptIndex[k]) {
            // there is a hole between (k-1) and k.
            wMuDt[ptIndex[k - 1] + 1] = (t[ptIndex[k - 1] + 1] - t[ptIndex[k - 1]]) /
                                        (t[ptIndex[k]] - t[ptIndex[k - 1]]);
            tem = (t[ptIndex[k - 1] + 1] - t[ptIndex[k - 1]]) *
                  (t[ptIndex[k]] - t[ptIndex[k - 1] + 1])
                  / (t[ptIndex[k]] - t[ptIndex[k - 1]]);

            // Test for NaN (z != z); 0/0 gives a NaN
            if (tem < 0 || tem != tem) {
               System.out.printf ("t[ptIndex[k-1]+1] - t[ptIndex[k-1]] = %g%n", t[ptIndex[k - 1] + 1] - t[ptIndex[k - 1]]);
               System.out.printf ("t[ptIndex[k]] - t[ptIndex[k-1]+1] = %g%n", t[ptIndex[k]] - t[ptIndex[k - 1] + 1]);
               System.out.printf ("t[ptIndex[k]] - t[ptIndex[k-1]] = %g%n", t[ptIndex[k]] - t[ptIndex[k - 1]]);
               System.out.printf ("t[ptIndex[k]] = %20.16g%n", t[ptIndex[k]]);
               System.out.printf ("t[ptIndex[k-1]] = %20.16g%n", t[ptIndex[k - 1]]);
               throw new IllegalStateException("   tem < 0 or NaN");
            }
            wSqrtDt[ptIndex[k - 1] + 1] = sigma * Math.sqrt (tem);
            wIndexList[indexCounter] = ptIndex[k] - 2;
            wIndexList[indexCounter + 1] = ptIndex[k] - 1;
            wIndexList[indexCounter + 2] = ptIndex[k];
            indexCounter += 3;
         }
      }
   }
}
